Análise Acústica da Igreja Matriz de Nossa Senhora da Conceição - Sabará (06/03/2017)

Preliminares

O bloco a seguir carrega as dependências necessárias. Estão são:

  • numpy - cálculo numérico e álgebra linear
  • scipy.fftpack - análise de Fourier
  • scipy.signal - processamento de sinais
  • scipy.stats - operações estatísticas
  • scipy.io.wavefile - lê e escreve arquivos .wav
  • matplotlib - gráficos
  • IPython - apenas para escutar os arquivos de áudio gerados no notebook
In [1]:
%matplotlib inline 
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.ticker
import IPython
import scipy.signal
import scipy.stats
from scipy.fftpack import fft, ifft, rfft
from scipy.io.wavfile import read as audioread
from scipy.io.wavfile import write as audiowrite
import seaborn as sns

E abaixo a função de geração do TSP do Hani:

In [2]:
from gen_tsp2 import gen_tsp2

Uso:

tsp, t = gen_tsp2(T=1200, fs=44.1, bw=22.05, bs=0, ta=120, tb=50, plot=False)

Parâmetros:

Parâmetro Descrição Unidade Valor Padrão
T Duração do TSP ms 1200
fs Taxa de amostragem kHz 44.1
bw Largura da faixa de passagem kHz 22.05
bs Deslocamento da faixa de passagem kHz 0
ta Início do TSP s 120
tb Taxa de crescimento do atraso de grupo ms/kHz 50
plot Mostrar os gráficos do TSP obtido - False

Valores de retorno:

Valor Descrição
tsp TSP no domínio do tempo
t Vetor temporal associado a x

Geração do TSP e obtenção da resposta em frequência

O bloco a seguir gera um TSP de 0 a 22,05 kHz com duração de 5s e o salva como um arquivo .wav:

In [3]:
# Parâmetros do TSP
T = 5000
fs = 44.1
bw = fs / 2
bs = 0
ta = T / 10
tb = 200

# Atribui o TSP gerado à variável 'tsp' e o vetor temporal associado a 't'
# Troque para plot=True se quiser ver o gráfico do TSP
tsp, t = gen_tsp2(T, fs, bw, bs, ta, tb, plot=False)

t /= 1000

# Salva o TSP gerado como um arquivo .wav
audiowrite(filename='tsp_python.wav',
           rate=int(fs*1e3),
           data=tsp.astype(np.float32))

Após reproduzir o TSP e gravar a resposta da sala, os arquivos produzidos são carregados:

In [4]:
# A lista abaixo contém os nomes dos arquivos a serem carregados, entre
# aspas e separados por vírgulas
lista_arquivos = ['exportadas/P1Mic1g80 - 2017-06-03 - Nossa Senhora da Conceição - 7s.wav',
                  'exportadas/P1Mic2g80 - 2017-06-03 - Nossa Senhora da Conceição - 7s.wav',
                  #'exportadas/P1Mic3g80 - 2017-06-03 - Nossa Senhora da Conceição - 7s.wav',
                  #'exportadas/P1Mic4g80 - 2017-06-03 - Nossa Senhora da Conceição - 7s.wav',
                  #'exportadas/P1Mic5g80 - 2017-06-03 - Nossa Senhora da Conceição - 7s.wav',
                  #'exportadas/P1Mic6g80 - 2017-06-03 - Nossa Senhora da Conceição - 7s.wav',
                 ]

# Após o bloco "for" a seguir, a variável 'respostas' conterá uma lista
# das respostas ao tsp carregadas
respostas = []
for arquivo in lista_arquivos:
    rate, resp = audioread(arquivo) # carrega a resposta do arquivo .wav
    resp = resp / np.max(np.abs(resp)) * .98 # normaliza para estar entre -1 e 1
    respostas.append(resp[:len(t)]) # põe na lista de respostas

# num_respostas é o número de respostas carregadas
num_respostas = len(respostas)

# players de áudio para ouvir as respostas carregadas
audio0 = IPython.display.Audio(data=tsp, rate=(fs * 1e3)) # tsp
label = IPython.display.Markdown('TSP')
IPython.display.display(label, audio0)

for indice, resposta in enumerate(respostas):
    titulo = IPython.display.Markdown(data='Resposta {}'.format(indice + 1)) # título
    audio = IPython.display.Audio(data=resposta, rate=(fs * 1e3)) # áudio
    IPython.display.display(titulo)
    IPython.display.display(audio)

/usr/local/lib/python3.6/site-packages/scipy/io/wavfile.py:273: WavFileWarning: Chunk (non-data) not understood, skipping it.
  WavFileWarning)

TSP

Resposta 1

Resposta 2

As respostas obtidas no domínio do tempo estão ilustradas abaixo:

In [5]:
# cria gráficos suficientes para plotar o TSP e todas as respostas
fig, graficos = plt.subplots(num_respostas + 1, figsize=(15, 20)) # 'graficos' é uma lista de gráficos
plt.tight_layout() # apenas estético

# plota o TSP no primeiro gráfico (0)
graficos[0].plot(t, tsp) # plota o tempo no eixo horizontal e tsp no vertical
graficos[0].set_title('TSP') # título do gráfico do TSP

# plota cada resposta nos graficos seguintes
for indice, resposta in enumerate(respostas): # para cada resposta
    graficos[indice + 1].plot(t, resposta) # plota o gráfico da resposta
    graficos[indice + 1].set_title('Resposta {}'.format(indice + 1)) # título da resposta

# insere o texto dos eixos
fig.text(.5, .001, 'Time (ms)', ha='center')
fig.text(-.02, .5, 'Normalized Amplitude', va='center', rotation='vertical')
#fig.savefig('Respostas_no_tempo.pdf', format='pdf')
Out[5]:
<matplotlib.text.Text at 0x10a846470>

Abaixo é feita a FFT de cada resposta:

In [6]:
# calcula um vetor de frequências 'f' (análogo ao vetor 't')
f = np.arange(0, len(t), 1) * (fs*1000) / len(t)

# calcula a FFT das respostas e as coloca em uma lista
fft_respostas = [fft(resp) for resp in respostas]

# e as respostas em frequência resultantes são passadas pra decibel (20 x log10 (abs(x))
db_fft_respostas = [20 * np.log10(np.abs(resp)) for resp in fft_respostas]

As respostas em frequência obtidas estão plotadas abaixo:

In [7]:
# cria uma lista das bandas padrão de 1/3 de oitava em kHz
freqs = np.array([.020, .025, .032, .040, .050, .063, .080, .100, .125, .160, .200, .250, .315, .400, .500, .630, .800,
                  1.0, 1.25, 1.6, 2.0, 2.5, 3.15, 4.0, 5.0, 6.3, 8.0, 10.0, 12.5, 16.0])

# cria gráficos suficientes para plotar o TSP e todas as respostas
fig, graficos = plt.subplots(num_respostas + 1, figsize=(15, 20)) # 'graficos' é uma lista de gráficos
plt.tight_layout() # apenas estético

# calcula a resposta em frequência do TSP para plotar (deve ser uma reta)
fft_tsp = fft(tsp)
db_fft_tsp = 20 * np.log10(np.abs(fft_tsp))

# plota o TSP no primeiro gráfico (0)
graficos[0].plot(f, db_fft_tsp) # plota as frequências no eixo horizontal e o tsp no vertical
graficos[0].set_title('TSP') # título do gráfico do TSP
graficos[0].set_xscale('log', basex=2)
graficos[0].set_xticks(1000 * freqs)
graficos[0].xaxis.set_major_formatter(matplotlib.ticker.ScalarFormatter())
graficos[0].axis([20, (fs*1000)/2, 0, 60]) # ajusta os limites do gráfico

# plota cada resposta nos graficos seguintes
for indice, resposta in enumerate(db_fft_respostas): # para cada resposta em frequência
    graficos[indice + 1].plot(f, resposta) # plota o gráfico da resposta
    graficos[indice + 1].set_title('Resposta {}'.format(indice + 1)) # título da resposta
    graficos[indice + 1].axis([20, (fs*1000)/2, 0, 60]) # ajusta os limites dos eixos
    
    # escala logarítmica
    graficos[indice + 1].set_xscale('log', basex=2)
    graficos[indice + 1].set_xticks(1000 * freqs)
    graficos[indice + 1].xaxis.set_major_formatter(matplotlib.ticker.ScalarFormatter())
    
    # põe linhas em todos os dós e fás
#     for i in range(0, 10):
#         c_freq = 32.7032 * 2**i
#         f_freq = 43.6535 * 2**i
#         graficos[indice + 1].axvline(c_freq, color='red', ls='-', label='Dó')
#         graficos[indice + 1].axvline(f_freq, color='black', ls='-', label='Fá')
        

# insere o texto dos eixos
fig.text(.5, .001, 'Frequency (Hz)', ha='center')
fig.text(-.02, .5, 'Amplitude (dB)', va='center', rotation='vertical')
#fig.savefig('Resposta_em_frequencia.pdf', format='pdf')
Out[7]:
<matplotlib.text.Text at 0x110e10940>

Derivação da resposta ao impulso a partir da resposta em frequência

A função 'get_ir_from_fr' abaixo calcula a resposta ao impulso da sala a partir de uma resposta em frequência. O cálculo é feito no domínio da frequência, i.e., se a resposta da sala que se deseja obter é $H$, o TSP é $X$ e a resposta gravada é $Y$, pode-se obter $H$ por $$H=\frac{Y}{X}=YX^{-1}$$

In [8]:
def get_ir_from_fft(fft_resp, fft_tsp):
    fft_ir = fft_resp / fft_tsp # deconvolui o TSP do sinal gravado
    ir = np.real(ifft(fft_ir)) # obtém a resposta ao impulso através da transformada inversa
    ir = ir / np.max(np.abs(ir)) * 0.98 # normalização
    
    return ir

Uso:

ir = get_ir_from_fft(fft_resp, fft_tsp)

Parâmetros:

Parâmetro Descrição
fft_resp FFT da resposta da sala ao TSP
fft_tsp FFT do TSP utilizado

Valores de retorno:

Valor Descrição
ir Vetor com a resposta ao impulso calculada

A linha a seguir calcula as respostas ao impulso a partir das respostas em frequência:

In [9]:
irs = [get_ir_from_fft(fft_resp, fft_tsp) for fft_resp in fft_respostas]

As respostas ao impulso obtidas estão plotadas abaixo:

In [10]:
# cria gráficos suficientes para plotar o TSP e todas as respostas
fig, graficos = plt.subplots(num_respostas + 1, figsize=(15, 20)) # 'graficos' é uma lista de gráficos
plt.tight_layout() # apenas estético

# plota o TSP no primeiro gráfico (0)
graficos[0].plot(t, tsp) # plota as frequências no eixo horizontal e o tsp no vertical
graficos[0].set_title('TSP') # título do gráfico do TSP

# plota cada resposta nos graficos seguintes
for indice, resposta in enumerate(irs): # para cada resposta ao impulso
    graficos[indice + 1].plot(t, resposta) # plota o gráfico da resposta
    graficos[indice + 1].set_title('Resposta {}'.format(indice + 1)) # título da resposta
    graficos[indice + 1].axis([0, 6, -1, 1])

# insere o texto dos eixos
fig.text(.5, .01, 'Time (ms)', ha='center')
fig.text(-.02, .5, 'Normalized Amplitude', va='center', rotation='vertical')
#fig.savefig('Resposta_ao_impulso.pdf', format='pdf')

# players de áudio para ouvir as respostas calculadas
for indice, resposta in enumerate(irs):
    titulo = IPython.display.Markdown(data='Resposta {}'.format(indice + 1)) # título
    audio = IPython.display.Audio(data=resposta, rate=(fs * 1e3)) # áudio
    IPython.display.display(titulo)
    IPython.display.display(audio)

Resposta 1

Resposta 2

Cálculo do Tempo de Reverberação (RT60)

A partir das respostas ao impulso obtidas, é possível calcular o tempo de reverberação para cada faixa de frequências. O cálculo é feito da seguinte forma:

  1. Para cada frequência, filtra-se o sinal. Foi utilizado para isso um filtro passa-faixas de butterworth de terceira ordem, com uma faixa de um terço de oitava centrado na frequência de interesse.
  2. Realiza-se a integração de Schroeder do sinal de forma a reduzir as variações, como descrito no ISO 8832-1. A integral é calculada da seguinte forma: $$S(t) = \int_{\infty}^{t}h^2(\tau)(-d\tau).$$ Para reduzir o impacto do ruído, o ISO 3382-1 recomenda realizar a integração até o ponto em que a curva decai ao nível do ruído de fundo. Assim, a operação a ser realizada é: $$S(t) = \int^{t}_{t1}h^2(\tau)(-d\tau),$$ onde t1 é o limite estimado. Este ponto foi estimado calculando a média móvel do sinal e estabelecendo o nível médio do último décimo do sinal como ruído.
  3. É feita uma regressão linear para se estimar o tempo que o sinal integrado leva para cruzar os pontos de -5 dB a -25 dB. Este tempo é o RT20.
  4. A partir do RT20, pode ser obter o RT60 como $$RT_{60}=3RT_{20}$$ As operações descritas são realizadas na função 'rt60', definida a seguir:
In [11]:
def rt60(ir, freqs, fs, beg=-5, end=-25, mul=3, plot=False):
    rt60 = [] # lista que conterá o RT60 para cada frequência pedida
    
    if plot:
        fig, axs = plt.subplots(len(freqs))
        if len(freqs) == 1:
            axs = [axs]
    
    for index, freq in enumerate(freqs): # para cada frequência
        # calcula as frequências de corte do filtro
        freq *= 1000 # frequência em Hz
        low = freq / 2**(1/6) # 2**(1/6) -> 1 terço de oitava
        low = low / (fs/2 * 1e3) # normalização
        high = freq * 2**(1/6)
        high = high / (fs/2 * 1e3)
            
        # calcula os parâmetros do filtro
        num, den = scipy.signal.butter(N=3, Wn=[low, high], btype='bandpass') # passa-faixas de terceira ordem
        
        # filtra o sinal e o normaliza
        filtered = scipy.signal.lfilter(num, den, ir)
        filtered = np.abs(filtered) / np.max(np.abs(filtered))
        
        # obtém a representação do sinal em dB para encontrar o limite de integração
        db_signal = 20 * np.log10(filtered)
        # obtém o início do impulso
        impulse = np.where(filtered == np.max(filtered))[0][0]
        # encontra o ponto em que o sinal atravessa o nível de ruído
        window = int(fs * 50) # 50 ms
        avg_signal = np.convolve(db_signal, np.ones(window) / window, mode='valid') # média móvel
        # estende o vetor da média móvel
        padding = int((len(db_signal) - len(avg_signal)) / 2)
        avg_signal = np.pad(avg_signal, padding, 'edge')
        threshold = int(len(db_signal) * 0.8) # último décimo do sinal
        noise_level = np.average(avg_signal[threshold:]) + 5 # calcula o nível de ruído como a média dos 10% finais
        
        t1 = np.abs(avg_signal[impulse:] - noise_level).argmin() + impulse # cruzamento do nível de ruído
        
        # realiza a integração de Schroeder
        sch = np.cumsum(filtered[t1:0:-1] ** 2)[::-1] # sinal integrado
        sch_db = 10 * np.log10(sch / np.max(sch)) # sinal integrado em dB
        
        # Encontra os pontos de cruzamento para a regressão
        #beg = -5 # inicio em -5 dB (ISO 3382-1)
        #end = -25 # -25 dB para RT20
        #mul = 3 # 3*RT20 = RT60
        first = sch_db[np.abs(sch_db - beg).argmin()] # cruzamento com -5 dB
        first_sample = np.where(sch_db == first)[0][0] # sample em que o cruzamento ocorre
        last = sch_db[np.abs(sch_db - end).argmin()] # cruzamento com -25 dB
        last_sample = np.where(sch_db == last)[0][0]
        
        # realiza a regressão linear nos pontos encontrados
        time = np.arange(first_sample, last_sample + 1) / fs
        sig = sch_db[first_sample:last_sample + 1]
        slope, intercept, r_value, p_value, std_err = scipy.stats.linregress(time, sig)
        reg_init = (beg - intercept) / slope # sample em que a regressão cruza -5 dB
        reg_end = (end - intercept) / slope # sample em que a regressão cruza -25 dB
        rt60.append(mul * (reg_end - reg_init))
        
        if plot:
            line = slope*np.arange(0, last_sample / fs + 1000) + intercept
            t = np.arange(0, len(db_signal)) / fs
            axs[index].plot(line)
            axs[index].plot(t, db_signal)
            axs[index].axvline(x=reg_init)
            axs[index].axvline(x=reg_end)
            axs[index].axhline(y=beg)
            axs[index].axhline(y=end)
            plt.show()
        
    return rt60

Uso:

lista_rt60 = rt60(ir, freqs, fs, plot=False)

Parâmetros:

Parâmetro Descrição Unidade Valor Padrão
ir Resposta ao impulso da sala - -
freqs Lista com as frequências desejadas kHz -
fs Taxa de amostragem kHz -
plot Exibir gráfico para cada frequência - False

Valores de retorno:

Valor Descrição Unidade
lista_rt60 Lista com o RT60 para cada frequência ms

O bloco a seguir calcula o RT60 as bandas padrão de 1/3 de oitava e os exibe em um gráfico, para cada resposta ao impulso obtida.

In [12]:
# cria uma lista das bandas padrão de 1/3 de oitava em kHz
freqs = np.array([.063, .080, .100, .125, .160, .200, .250, .315, .400, .500, .630, .800,
                  1.0, 1.25, 1.6, 2.0, 2.5, 3.15, 4.0, 5.0, 6.3, 8.0, 10.0, 12.5, 16.0])

# calcula e plota os RT60 para cada frequência
lista_edt = list()
lista_rt60_por_rt20 = list() # rt60 calculado para cada amostra
lista_rt60_por_rt30 = list()
lista_rt60_por_rt60 = list()
for resposta in irs: # para cada resposta ao impulso
    lista_edt.append(np.array(rt60(resposta, freqs, fs, beg=0, end=-10, mul=1))) # obtém o edt
    lista_rt60_por_rt20.append(np.array(rt60(resposta, freqs, fs, beg=-5, end=-25, mul=3))) # obtém o RT60 pelo RT20
    lista_rt60_por_rt30.append(np.array(rt60(resposta, freqs, fs, beg=-5, end=-35, mul=2))) # obtém o RT60 pelo RT30
    lista_rt60_por_rt60.append(np.array(rt60(resposta, freqs, fs, beg=-5, end=-65, mul=1))) # obtém o RT60 diretamente
In [13]:
# cria uma lista das bandas padrão de 1/3 de oitava em kHz
freqs = np.array([.063, .080, .100, .125, .160, .200, .250, .315, .400, .500, .630, .800,
                  1.0, 1.25, 1.6, 2.0, 2.5, 3.15, 4.0, 5.0, 6.3, 8.0, 10.0, 12.5, 16.0])

# cria um número suficiente de gráficos
fig, axs = plt.subplots(num_respostas + 1, figsize=(15, 20)) # um gráfico adicional para a média
plt.tight_layout()

if num_respostas == 1:
    axs = [axs]

# calcula e plota os RT60 para cada frequência
for indice in range(len(irs)): # para cada resposta ao impulso
    axs[indice].plot(freqs*1000, np.array([
                                           #lista_edt[indice] / 1000,
                                           lista_rt60_por_rt20[indice] / 1000, 
                                           #lista_rt60_por_rt30[indice] / 1000, 
                                           #lista_rt60_por_rt60[indice] / 1000,
                                          ]).transpose()) # plota os tempos obtidos em um gráfico
    axs[indice].set_title('Resposta {}'.format(indice + 1))
    axs[indice].axis([64, fs*1000/2, 0, 7]) # alterar os limites de acordo com a situação
    
    # linha mediana
    mediana = np.median(lista_rt60_por_rt20[indice])
    rt60_1k = lista_rt60_por_rt20[indice][12] /1000
    axs[indice].axhline(y=rt60_1k, color='red', ls='--', label='RT60 em 1.000 Hz ={:.2f} s'.format(rt60_1k))
    
    # escala logarítmica
    axs[indice].set_xscale('log', basex=2)
    axs[indice].set_xticks(1000 * freqs)
    axs[indice].xaxis.set_major_formatter(matplotlib.ticker.ScalarFormatter())
    
    # legenda
    handles, labels = axs[indice].get_legend_handles_labels()
    axs[indice].legend(handles, labels)
    
# calcula a média e o desvio padrão das amostras em cada frequência
media_rt60, desvio_rt60 = zip(*[(freq.mean() / 1000, freq.std() / 1000) for freq in np.array(lista_rt60_por_rt20).transpose()])
#axs[-1].plot(freqs*1000, np.array(media_rt60))
axs[-1].errorbar(freqs*1000, media_rt60, yerr=desvio_rt60, ecolor='black', label='RT60 com desvio padrão')
axs[-1].set_title('RT60 médio da igreja')
axs[-1].axis([64, fs*1000/2, 0, 2.5]) # alterar os limites de acordo com a situação

# linha mediana
mediana = np.median(media_rt60) / 1000
rt60_1k = media_rt60[12]
axs[-1].axhline(y=rt60_1k, color='red', ls='--', label='RT60 em 1.000 Hz = {:.2f} s'.format(rt60_1k))

# escala logarítmica
axs[-1].set_xscale('log', basex=2)
axs[-1].set_xticks(1000 * freqs)
axs[-1].xaxis.set_major_formatter(matplotlib.ticker.ScalarFormatter())

# legenda
handles, labels = axs[-1].get_legend_handles_labels()
axs[-1].legend(handles, labels)
    
# insere o texto dos eixos
fig.text(.5, .001, 'Frequency (Hz)', ha='center')
fig.text(-.02, .5, 'RT60 (s)', va='center', rotation='vertical')
#fig.savefig('RT60.pdf', format='pdf')
Out[13]:
<matplotlib.text.Text at 0x1124fbfd0>

Cálculo dos demais parâmetros acústicos

Claridade (C50 e C80)

Os parâmetros de claridade são a razão entre a energia sonora antes e depois de um determinado instante, i.e. $$C_t=10 \log_{10}\left( \frac{ \int_0^t p^2(\tau) d\tau }{ \int_t^\infty p^2(\tau) d\tau } \right).$$ Os valores mais utilizados são 50 ms (correlacionado com a inteligibilidade da fala) e 80 ms (para música). A função abaixo calcula a claridade para qualquer valor de tempo e determinadas frequências, a partir da resposta ao impulso:

In [14]:
def clarity(ir, time, freqs, fs):
    init_sample = np.where(ir == np.max(ir))[0][0] # sample do início do impulso
    c = [] # lista onde serão armazenados os valores de claridade
    
    for freq in freqs: # para cada frequência
        # obtém as frequências de corte do filtro
        freq *= 1e3 # frequencia em Hz
        low = freq / 2**(1/6) # 1/3 de oitava
        low = low / (fs/2 * 1e3) # normalização
        high = freq * 2**(1/6)
        high = high / (fs/2 * 1e3)
        
        # calcula os parâmetros do filtro
        num, den = scipy.signal.butter(3, [low, high], 'bandpass')
        
        # filtra o sinal na faixa desejada
        filtered_signal = scipy.signal.lfilter(num, den, ir)
        
        # calcula a claridade pela fórmula acima
        h2 = filtered_signal ** 2 # quadrado do sinal
        t = init_sample + int(time*fs + 1) # valor de t na fórmula
        c.append(10 * np.log10(np.sum(h2[init_sample:t]) / np.sum(h2[t:]))) # integração
        
    return c

Uso:

lista_claridades = clarity(ir, time, freqs, fs)

Parâmetros:

Parâmetro Descrição Unidade
ir Resposta ao impulso calculada -
time Instante desejado, ex. 50 para C50 ms
freqs Lista de frequências kHz
fs Taxa de amostragem kHz

Valor retornado:

Valor Descrição
lista_claridades Lista com o valor da claridade para o tempo escolhido em cada frequência

O bloco abaixo calcula os valores de C50 para várias frequências das respostas obtidas e os plota em gráficos.

In [15]:
# cria uma lista das bandas padrão de 1/3 de oitava em kHz
freqs = np.array([.063, .080, .100, .125, .160, .200, .250, .315, .400, .500, .630, .800,
                  1.0, 1.25, 1.6, 2.0, 2.5, 3.15, 4.0, 5.0, 6.3, 8.0, 10.0, 12.5, 16.0])

# cria um número suficiente de gráficos
fig, axs = plt.subplots(num_respostas + 1, figsize=(15, 20))
plt.tight_layout()

if num_respostas == 1:
    axs = [axs]

# plota os C50 para cada frequência
lista_c50 = list()
for indice, resposta in enumerate(irs): # para cada resposta ao impulso
    lista_c50.append(clarity(resposta, 50, freqs, fs)) # obtém o C50
    axs[indice].plot(freqs*1000, lista_c50[indice]) # plota os tempos obtidos em um gráfico
    axs[indice].set_title('Resposta {}'.format(indice + 1))
    axs[indice].axis([min(freqs)*1000, fs*1000/2, -10, 15]) # alterar os limites de acordo com a situação
    
    # linha mediana
    mediana = np.median(lista_c50[indice])
    c50_1k = lista_c50[indice][12]
    axs[indice].axhline(y=c50_1k, color='red', ls='--', label='C50={:.2f} dB'.format(c50_1k))
    
    # escala logarítmica
    axs[indice].set_xscale('log', basex=2)
    axs[indice].set_xticks(1000 * freqs)
    axs[indice].xaxis.set_major_formatter(matplotlib.ticker.ScalarFormatter())
    
    # legenda
    handles, labels = axs[indice].get_legend_handles_labels()
    axs[indice].legend(handles, labels)
    
# calcula a média e o desvio padrão das amostras em cada frequência
media_c50, desvio_c50 = zip(*[(freq.mean(), freq.std()) for freq in np.array(lista_c50).transpose()])
#axs[-1].plot(freqs*1000, media_c50)
axs[-1].errorbar(freqs*1000, media_c50, yerr=desvio_c50)
axs[-1].set_title('Média')
# escala logarítimica
axs[-1].set_xscale('log', basex=2)
axs[-1].set_xticks(1000 * freqs)
axs[-1].xaxis.set_major_formatter(matplotlib.ticker.ScalarFormatter())
# linha mediana
mediana = np.median(media_c50)
c50_1k = media_c50[12]
axs[-1].axhline(y=c50_1k, color='red', ls='--', label='C50={:.2f} dB'.format(c50_1k))
axs[-1].axis([min(freqs) * 1000, fs*1000/2, -10, 10]) # limites
# legenda
handles, labels = axs[-1].get_legend_handles_labels()
axs[-1].legend(handles, labels)

# insere o texto dos eixos
fig.text(.5, .001, 'Frequency (Hz)', ha='center')
fig.text(-.02, .5, 'C50 (dB)', va='center', rotation='vertical')
#fig.savefig('C50.pdf', format='pdf')
Out[15]:
<matplotlib.text.Text at 0x113ae9da0>

E o bloco abaixo faz o mesmo com o C80 de cada resposta

In [16]:
# cria uma lista das bandas padrão de 1/3 de oitava em kHz
freqs = np.array([.063, .080, .100, .125, .160, .200, .250, .315, .400, .500, .630, .800,
                  1.0, 1.25, 1.6, 2.0, 2.5, 3.15, 4.0, 5.0, 6.3, 8.0, 10.0, 12.5, 16.0])

# cria um número suficiente de gráficos
fig, axs = plt.subplots(num_respostas + 1, figsize=(15, 20))
plt.tight_layout()

if num_respostas == 1:
    axs = [axs]

# plota os C80 para cada frequência
lista_c80 = list()
for indice, resposta in enumerate(irs): # para cada resposta ao impulso
    lista_c80.append(clarity(resposta, 80, freqs, fs)) # obtém o C80
    axs[indice].plot(freqs*1000, lista_c80[indice]) # plota os tempos obtidos em um gráfico
    axs[indice].set_title('Resposta {}'.format(indice + 1))
    axs[indice].axis([min(freqs)*1000, fs*1000/2, -17.5, 15]) # alterar os limites de acordo com a situação
    axs[indice].fill_between(freqs*1000, -2, 2, color='green', alpha=.25,label='Faixa de variacão ótima (-2dB à +2dB)')
    axs[indice].fill_between(freqs*1000, -5, -2, color='green', alpha=.1, label='Faixa de variacão típica (-5dB à +5dB)')
    axs[indice].fill_between(freqs*1000, 2, 5, color='green', alpha=.1)
    
    # linha mediana
    mediana = np.median(lista_c80[indice])
    c80_1k = lista_c80[indice][12]
    axs[indice].axhline(y=c80_1k, color='red', ls='--', label='C80 em 1.000Hz = {:.2f} dB'.format(c80_1k))
    
    
    # escala logarítmica
    axs[indice].set_xscale('log', basex=2)
    axs[indice].set_xticks(1000 * freqs)
    axs[indice].xaxis.set_major_formatter(matplotlib.ticker.ScalarFormatter())
    
    # legenda
    handles, labels = axs[indice].get_legend_handles_labels()
    axs[indice].legend(handles, labels, loc='lower right')
    
# calcula a média e o desvio padrão das amostras em cada frequência
media_c80, desvio_c80 = zip(*[(freq.mean(), freq.std()) for freq in np.array(lista_c80).transpose()])
#axs[-1].plot(freqs*1000, media_c80)
axs[-1].errorbar(freqs*1000, media_c80, yerr=desvio_c80, label='C80 com desvio padrão')
axs[-1].set_title('C80 médio da igreja')
axs[-1].fill_between(freqs*1000, -2, 2, color='green', alpha=.25, label='Faixa de variacão ótima (-2dB à +2dB)')
axs[-1].fill_between(freqs*1000, -5, -2, color='green', alpha=.1, label='Faixa de variacão típica (-5dB à +5dB)')
axs[-1].fill_between(freqs*1000, 2, 5, color='green', alpha=.1)

# escala logarítimica
axs[-1].set_xscale('log', basex=2)
axs[-1].set_xticks(1000 * freqs)
axs[-1].xaxis.set_major_formatter(matplotlib.ticker.ScalarFormatter())

# linha mediana
mediana = np.median(media_c80)
c80_1k = media_c80[12]
axs[-1].axhline(y=c80_1k, color='red', ls='--', label='C80 em 1.000 Hz = {:.2f} dB'.format(c80_1k))
axs[-1].axis([min(freqs) * 1000, fs*1000/2, -17.5, 15]) # limites
    

# legenda
handles, labels = axs[-1].get_legend_handles_labels()
axs[-1].legend(handles, labels, loc='lower right')

# insere o texto dos eixos
fig.text(.5, .001, 'Frequency (Hz)', ha='center')
fig.text(-.02, .5, 'C80 (dB)', va='center', rotation='vertical')
#fig.savefig('C80.pdf', format='pdf')
Out[16]:
<matplotlib.text.Text at 0x114603278>

Definição (D50)

A definição é um outro indicador da inteligibilidade da fala. É a razão entre a energia sonora nos primeiros 50 ms e a de todo o período de captura, i.e. $$D_{50} = \frac{ \int_0^{50} p^2(\tau) d\tau }{ \int_0^\infty p^2(\tau) d\tau }.$$ A função abaixo calcula a claridade de uma resposta ao impulso em várias frequências.

In [17]:
def definition(ir, freqs, fs):
    init_sample = np.where(ir == np.max(ir))[0][0] # amostra do início do impulso
    d = [] # lista onde serão armazenados os valores de definição
    
    for freq in freqs: # para cada frequência
        # obtém as frequências de corte do filtro
        freq *= 1e3 # frequencia em Hz
        low = freq / 2**(1/6) # 1/3 de oitava
        low = low / (fs/2 * 1e3) # normalização
        high = freq * 2**(1/6)
        high = high / (fs/2 * 1e3)
        
        # calcula os parâmetros do filtro
        num, den = scipy.signal.butter(3, [low, high], 'bandpass')
        
        # filtra o sinal na faixa desejada
        filtered_signal = scipy.signal.lfilter(num, den, ir)
        
        # calcula o D50 pela fórmula acima
        h2 = filtered_signal ** 2 # quadrado do sinal
        t = init_sample + int(50*fs + 1) # 50 ms após o início do impulso
        d.append(np.sum(h2[init_sample:t]) / np.sum(h2[init_sample:])) # integração
        
    return d

Uso:

lista_definicoes = definition(ir, freqs, fs)

Parâmetros:

Parâmetro Descrição Unidade
ir Resposta ao impulso calculada -
freqs Lista de frequências kHz
fs Taxa de amostragem kHz

Valor retornado:

Valor Descrição
lista_definitions Lista com o valor da definição da resposta ir em cada frequência

O bloco a seguir calcula a definição em várias frequências de cada uma das respostas obtidas e plota o resultado em um gráfico.

In [ ]:
In [18]:
# cria uma lista das bandas padrão de 1/3 de oitava em kHz
freqs = np.array([.063, .080, .100, .125, .160, .200, .250, .315, .400, .500, .630, .800,
                  1.0, 1.25, 1.6, 2.0, 2.5, 3.15, 4.0, 5.0, 6.3, 8.0, 10.0, 12.5, 16.0])                  

# cria um número suficiente de gráficos
fig, axs = plt.subplots(num_respostas + 1, figsize=(15, 20))
plt.tight_layout()

if num_respostas == 1:
    axs = [axs]

# plota os D50 para cada frequência
lista_d50 = list()
for indice, resposta in enumerate(irs): # para cada resposta ao impulso
    lista_d50.append(definition(resposta, freqs, fs)) # obtém o D50
    axs[indice].plot(freqs*1000, lista_d50[indice]) # plota os tempos obtidos em um gráfico
    axs[indice].set_title('Resposta {}'.format(indice + 1))
    axs[indice].axis([min(freqs)*1000, fs*1000/2, 0, 1]) # alterar os limites de acordo com a situação
    axs[indice].fill_between(freqs*1000, .3, .7, color='green', alpha=.2, label='Faixa de variação típica (0.3 à 0.7)')
    
    # linha mediana
    mediana = np.median(lista_d50[indice])
    d50_1k = lista_d50[indice][12]
    axs[indice].axhline(y=d50_1k, color='red', ls='--', label='D50={:.2f}'.format(d50_1k))
    
    # escala logarítmica
    axs[indice].set_xscale('log', basex=2)
    axs[indice].set_xticks(1000 * freqs)
    axs[indice].xaxis.set_major_formatter(matplotlib.ticker.ScalarFormatter())
    
    # legenda
    handles, labels = axs[indice].get_legend_handles_labels()
    axs[indice].legend(handles, labels, loc='lower right')
    
# calcula a média e o desvio padrão das amostras em cada frequência
media_d50, desvio_d50 = zip(*[(freq.mean(), freq.std()) for freq in np.array(lista_d50).transpose()])
#axs[-1].plot(freqs*1000, media_d50, ls='-')
#axs[-1].plot(freqs*1000, np.array(lista_d50).transpose()) # plota todas as respostas em um só gráfico
axs[-1].errorbar(freqs*1000, media_d50, yerr=desvio_d50, label='D50 com desvio padrão')
axs[-1].set_title('D50 médio da igreja')
axs[-1].fill_between(freqs*1000, .3, .7, color='green', alpha=.2, label='Faixa de variação típica (0.3 à 0.7)')

# escala logarítimica
axs[-1].set_xscale('log', basex=2)
axs[-1].set_xticks(1000 * freqs)
axs[-1].xaxis.set_major_formatter(matplotlib.ticker.ScalarFormatter())

# linha mediana
mediana = np.median(media_d50)
d50_1k = media_d50[12]
axs[-1].axhline(y=d50_1k, color='red', ls='--', label='D50 em 1.000 Hz = {:.2f}'.format(d50_1k))
axs[-1].axis([min(freqs) * 1000, fs*1000/2, 0, 1]) # limites

# legenda
handles, labels = axs[-1].get_legend_handles_labels()
axs[-1].legend(handles, labels,loc='lower right')

# insere o texto dos eixos
fig.text(.5, .001, 'Frequency (Hz)', ha='center')
fig.text(-.02, .5, 'D50', va='center', rotation='vertical')
#fig.savefig('D50.pdf', format='pdf')
Out[18]:
<matplotlib.text.Text at 0x130ce5780>

Centróide Espectral

O centróide espectral é uma indicação do "centro de massa" de um sinal. É uma média das frequências do sinal, ponderada por suas amplitudes, i.e., $$\text{Centroide} = \frac{\sum_{n=0}^{N-1}f(n)X(n)}{\sum_{n=0}^{N-1}X(n)}$$ onde f(n) é a n-ésima frequência central dada pela FFT e X(n) sua amplitude.

In [19]:
def centroid(ir, fs):
    X = np.abs(rfft(ir)) # magnitudes
    #X = X / sum(X) # normalização
    freqs = np.linspace(0, 1, len(ir))
    bins = np.abs(np.fft.fftfreq(len(ir), 1.0 / fs)) # frequencias da fft
    return np.sum(bins * X) / np.sum(X) #* bins[-1] # média ponderada
In [20]:
# cria uma lista das bandas padrão de 1/3 de oitava em kHz
freqs = np.array([.020, .025, .032, .040, .050, .063, .080, .100, .125, .160, .200, .250, .315, .400, .500, .630, .800,
                  1.0, 1.25, 1.6, 2.0, 2.5, 3.15, 4.0, 5.0, 6.3, 8.0, 10.0, 12.5, 16.0])

# cria gráficos suficientes para plotar todas as respostas
fig, graficos = plt.subplots(num_respostas + 1, figsize=(15,20)) # 'graficos' é uma lista de gráficos
plt.tight_layout() # apenas estético

lista_centroides = []

# plota cada resposta nos graficos seguintes
for indice, resposta in enumerate(db_fft_respostas): # para cada resposta em frequência
    graficos[indice].plot(f, resposta) # plota o gráfico da resposta
    graficos[indice].set_title('Resposta {}'.format(indice + 1)) # título da resposta
    graficos[indice].axis([16, (fs*1000)/2, 0, 60]) # ajusta os limites dos eixos
    
    # escala logarítmica
    axs[indice].set_xscale('log', basex=2)
    axs[indice].set_xticks(1000 * freqs)
    axs[indice].xaxis.set_major_formatter(matplotlib.ticker.ScalarFormatter())
    
    # centróide espectral
    centroide = centroid(resposta, fs*1000)
    lista_centroides.append(centroide)
    graficos[indice].axvline(centroide, color='red', ls='-', label='Centróide: {:.0f} Hz'.format(centroide))
    
    # legenda
    handles, labels = graficos[indice].get_legend_handles_labels()
    graficos[indice].legend(handles, labels)

# plota a média dos centróides
centroide_medio, desvio_centroide = np.mean(lista_centroides), np.std(lista_centroides)
resposta_media = sum(db_fft_respostas) / len(db_fft_respostas)
graficos[-1].plot(f, resposta_media)
graficos[-1].set_title('Centróide de frequência médio da igreja')
graficos[-1].axis([16, (fs*1000) / 2, 0, 60])
graficos[-1].axvline(centroide_medio, color='red', ls='-', label='Centróide: {:.0f} Hz'.format(centroide_medio))
handles, labels = graficos[-1].get_legend_handles_labels()
graficos[-1].legend(handles, labels)

# insere o texto dos eixos
fig.text(.5, .001, 'Frequency (Hz)', ha='center')
fig.text(-.02, .5, 'Amplitude (dB)', va='center', rotation='vertical')
#fig.savefig('centroide_de_frequencia.png', format='png')
Out[20]:
<matplotlib.text.Text at 0x14abafef0>

Auralização

In [21]:
rate, voice = audioread('Voice.16bit.wav')
_, musica_lenta = audioread('HarpMono.wav')
_, musica_rapida = audioread('musica_rapida.wav')

/usr/local/lib/python3.6/site-packages/scipy/io/wavfile.py:273: WavFileWarning: Chunk (non-data) not understood, skipping it.
  WavFileWarning)

Fala

In [22]:
auralizacao = np.convolve(voice, irs[1][:int(2e3*fs)])
In [23]:
titulo1 = IPython.display.Markdown(data='Original') # título
audio1 = IPython.display.Audio(data=voice, rate=rate) # áudio
IPython.display.display(titulo1)
IPython.display.display(audio1)
titulo = IPython.display.Markdown(data='Auralização') # título
audio = IPython.display.Audio(data=auralizacao, rate=rate) # áudio
IPython.display.display(titulo)
IPython.display.display(audio)
audiowrite(filename='auralizacao.wav',
           rate=rate,
           data=(0.9*auralizacao).astype(np.int16))

Original

Auralização

Música

In [25]:
auralizacao_musica_lenta = np.convolve(musica_lenta, irs[1][:int(2e3*fs)])
auralizacao_musica_rapida = np.convolve(musica_rapida, irs[1][:int(2e3*fs)])
In [26]:
titulo1 = IPython.display.Markdown(data='Música 1 - Original') # título
audio1 = IPython.display.Audio(data=musica_lenta, rate=rate) # áudio
IPython.display.display(titulo1)
IPython.display.display(audio1)
titulo2 = IPython.display.Markdown(data='Música 2 - Original') # título
audio2 = IPython.display.Audio(data=musica_rapida, rate=rate) # áudio
IPython.display.display(titulo2)
IPython.display.display(audio2)
titulo3 = IPython.display.Markdown(data='Música 1 - Auralização') # título
audio3 = IPython.display.Audio(data=auralizacao_musica_lenta, rate=rate) # áudio
IPython.display.display(titulo3)
IPython.display.display(audio3)
audiowrite(filename='auralizacao_musica_lenta.wav',
           rate=rate,
           data=(0.9*auralizacao_musica_lenta).astype(np.int16))
titulo4 = IPython.display.Markdown(data='Música 2 - Auralização') # título
audio4 = IPython.display.Audio(data=auralizacao_musica_rapida, rate=rate) # áudio
IPython.display.display(titulo4)
IPython.display.display(audio4)
audiowrite(filename='auralizacao_musica_rapida.wav',
           rate=rate,
           data=(0.9*auralizacao_musica_rapida).astype(np.int16))

Música 1 - Original

Música 2 - Original

Música 1 - Auralização

Música 2 - Auralização